from_loop = [1 1 2718 2718  1 ] ; % 1   
until_loop = [ 10000 6685 10000 6685 5020  ] ; %  

rounds = length(from_loop) ;
collect_base_irf = zeros(2,T_irf,length(choose_variab)) ;
collect_other_irf = zeros(length(choose_variab),T_irf,rounds) ;

for sample = 1:rounds


%% Organize DATA

names = fullnames_data(choose_variab)
bettername = betternames_data(choose_variab) 
region = region_data(choose_variab)
names_exog = fullnames_data(choose_exog)
Y = data_all(:,choose_variab);
Y_exog=data_all(:,choose_exog);

% chose names for IRFs
names_IRFs = fullnames_data(choose_variab);
transformIRF_used = transformIRFs_data(choose_variab);
labelIRF_used = labelIRFs_data(choose_variab);
Numbers = 1:length(day);
Numbers=Numbers';

% reduce to period for which all measures are available
cond = isnan(Y);
cond2 = sum(cond,2) == 0;
Y_start = Y(cond2,:);
Y_exog_start = Y_exog(cond2,:);
year_available     = year(cond2);
month_available    = month(cond2);
day_available    = day(cond2);
Numbers_available = Numbers(cond2);
% plot(Y(:,2));

data_all_available = data_all(cond2,:);
% keep only after a certain period, in case data are available for all variables added for before that period


from = from_loop(sample) ;
until = until_loop(sample) ;

bootstrap_on = and(from == 1, until == 10000 ); 

cond = and(Numbers_available >= from,Numbers_available <= until);
sum(cond)
Y = Y_start(cond,:);
Y_exog = Y_exog_start(cond,:);

year_available = year_available(cond);
month_available = month_available(cond);
day_available    = day_available(cond);
Numbers_available = Numbers_available(cond);

disp(['-->> Data available for :    ', num2str(year_available(1)) 'M' num2str(month_available(1)) 'D' num2str(day_available(1)) ...
                            ' until ', num2str(year_available(end)) 'M' num2str(month_available(end)) 'D' num2str(day_available(end))] );

% ensure Y enters horizontally
Y = Y';
Y_exog = Y_exog';
k = size(Y,1)
T = size(Y,2)
assert(  T > k)
assert( length(year_available) == T)      

year_available_help=zeros(T,1);
for i=1:length(year_available)
year_available_help(i,:) = year_available(i,:) + 0.001*i;
end


%% Estimate RVAR
maxlags = 20;
 
constantD   = 1;
trendD      = 0;
trendBreaks = [];
seasonD     = 0;
moreBreaks  = [];
 

[Ahat, Resid, Sigmahat, Ahat_determ, X_determ, T_used, df_lost] = ...
         zFC_VAR_Estimation_exog(Y,p,constantD,trendD,trendBreaks,seasonD,moreBreaks,exogD,Y_exog);        
        
year_available_used     = year_available(p+1:end);
assert(length(year_available_used) == size(Resid,2))
month_available_used    = month_available(p+1:end);
day_available_used    = day_available(p+1:end);
Numbers_available_used = Numbers_available(p+1:end);

%% IRFs, point estimation
T_step = [1:1:T_irf]; % used to plot IRFs
shock_considered = 1; % variable shocked

% CHOLESKY
if cholesky_on == 1
    B_chol          = chol(Sigmahat)';
    shock_size_chol = 1;
    IRF_chol        = zFC_IRFs(Ahat, B_chol, T_irf, shock_considered, shock_size_chol);
end 

Instrument = data_all(:,choose_shock);
cond_i = isnan(Instrument);
cond2_i = sum(cond_i,2) == 0;
Instrument = Instrument(cond2_i,:);
year_available_ins     = year(cond2_i);
month_available_ins    = month(cond2_i);
day_available_ins = day(cond2_i);
Numbers_available_ins = Numbers(cond2_i);

% select subperiod common to shocks used and variables used
from2_year         = max(year_available_used(1),year_available_ins(1));
from2_month          = max(month_available_used(1),month_available_ins(1));
from2_day          = max(day_available_used(1),day_available_ins(1));
from2          = max(Numbers_available_used(1),Numbers_available_ins(1));
% % % % %from2          = timeaxis_used(59);
until2_year         = min(year_available_used(end),year_available_ins(end));
until2_month         = min(month_available_used(end),month_available_ins(end));
until2_day         = min(day_available_used(end),day_available_ins(end));
until2         = min(Numbers_available_used(end),Numbers_available_ins(end));
Numers_available_used2 = [from2:until2];
%                         
disp(['-->> Period effectively used in the identification:    ', num2str(year_available_ins(1)) 'M' num2str(month_available_ins(1)) 'D' num2str(day_available_ins(1)) ...
                            ' until ', num2str(year_available_ins(end)) 'M' num2str(month_available_ins(end)) 'D' num2str(day_available_ins(end))] );
                      
% Take only the instruments in this subperiod
cond1            = Numbers_available_ins >= from2 & Numbers_available_ins <= until2;
Instruments_used = Instrument(cond1,:);
n_instrum=size(Instrument);
n_instrum=n_instrum(1,2);
 
% Take only the residuals in this subperiod
cond2      = Numbers_available_used >= from2 & Numbers_available_used <= until2;
Resid_used = Resid(:,cond2);
AAA=[Resid_used' Instruments_used];


names_instrument = fullnames_data(choose_shock) ;
[bla,g] = size(Instruments_used) ;

% Identify the model
% choose correct position for indicator variable
position_indicator = 1;

[imp_vect_relat, absolute_impulse_vector, ...
             F_stat_test1, p_values_test1, signif_test1, Rsquared_test1, Betas_test1, ...
             F_stat_test2, p_values_test2, signif_test2, Rsquared_test2, ...
             correl_instrument_shock, a_vect, estim_shocks] = zFC_IdentifExtInstruments5(Resid_used, Instruments_used, position_indicator, Sigmahat);
% 
table_validity_rows = char('beta', 'numb. of stars', 'numb. of obs.', 'F statistic', 'R squared', 'obs. of the proxy used');
table_validity_rows = cellstr(table_validity_rows); 
table_validity = [Betas_test1'; signif_test1(:,1)'; length(Instruments_used)*ones(1,k)/1000; F_stat_test1'; Rsquared_test1']    

nonzero_instrument = nonzeros(Instruments_used);
n_nonzero_instrument = length(nonzero_instrument)

B_extin    = [absolute_impulse_vector, zeros(k,k-1)];
shock_size_extin = shock_size/absolute_impulse_vector(1);
IRF_extin  = zFC_IRFs(Ahat, B_extin, T_irf, 1, shock_size_extin);

%% Bootstrap on the reduced form model 

if bootstrap_on == 1 
    
    Yinit           = Y(:,1:p);
    Ahat_loop       = NaN*ones(k,k*p,Loops);
    Resid_used_loop = NaN*ones(k,length(Resid_used),Loops);
    Sigmahat_loop   = NaN*ones(k,k,Loops);
    Instruments_used_loop = NaN*ones(length(Instruments_used),n_instrum,Loops);

    for l = 1:Loops

        % Generate a vector that flips signs randomly
        step = randn(T_used,1);
        change_sign = -1*ones(T_used,1);
        change_sign(step > 0) = 1;

        % Generate pseudo data
        Resid_pseudo = Resid*diag(change_sign);

        Y_pseudo = zFC_GenerateVARData(Yinit, Ahat, Resid_pseudo, Ahat_determ, X_determ);

        % Estimate the VAR on the pseudo data
        [Ahat_loop(:,:,l), Resid_step, Sigmahat_loop(:,:,l), ~, ~, ~, ~] = ...
                       zFC_VAR_Estimation_exog(Y_pseudo,p,constantD,trendD,trendBreaks,seasonD,moreBreaks,exogD,Y_exog);

        % Store residuals for the subpart that overlaps with the shocks used
        Resid_used_loop(:,:,l) = Resid_step(:,cond2);

        % Flip sign of the instruments using the same order in the bootstrapped residuals, consider only the part that overlaps with the residuals
        Instruments_used_loop(:,:,l) = diag(change_sign(cond2))*Instruments_used;
    end

    % IRFs, error bands
    IRF_chol_loop  = NaN*ones(k,T_irf,Loops);
    IRF_extin_loop = NaN*ones(k,T_irf,Loops);
    IRF_extin_new_loop = NaN*ones(k,T_irf,Loops);

    for l = 1:Loops

        % External instrument, using only monthly data
        [~, ~, relative_impulse_vector_loop, ~, incons_impulse_vector_loop] = zFC_IdentifExtInstruments2(Resid_used_loop(:,:,l),Instruments_used_loop(:,:,l));
        b11_loop = zFC_ScaleInExternalInstrument(Sigmahat_loop(:,:,l), relative_impulse_vector_loop);
        absolute_impulse_vector_loop = relative_impulse_vector_loop*b11_loop;

        B_extin_loop    = [absolute_impulse_vector_loop, zeros(k,k-1)];
         shock_size_extin_loop = shock_size/absolute_impulse_vector_loop(1);
        IRF_extin_loop(:,:,l)  = zFC_IRFs(Ahat_loop(:,:,l), B_extin_loop, T_irf, 1, shock_size_extin_loop);

    end

    % Compute error bands
    if cholesky_on ==1
        IRF_chol_bands = NaN*ones(k,T_irf,2);
    end

    IRF_extin_bands = NaN*ones(k,T_irf,2);
    interval       = 10;
    for ttt = 1:T_irf 
        for kkk = 1:k
            if cholesky_on == 1
                IRF_chol_bands(kkk,ttt,1)   =  prctile(IRF_chol_loop(kkk,ttt,:),interval/2);
                IRF_chol_bands(kkk,ttt,2)   =  prctile(IRF_chol_loop(kkk,ttt,:),100-interval/2);   
            end
            IRF_extin_bands(kkk,ttt,1)  =  prctile(IRF_extin_loop(kkk,ttt,:),interval/2);
            IRF_extin_bands(kkk,ttt,2)  =  prctile(IRF_extin_loop(kkk,ttt,:),100-interval/2);   
        end
    end

end

%% collect irfs

for i = 1:k
   
    collect_base_irf(1,:,i) = transformIRF_used(i)*IRF_extin_bands(i,:,2) ;
    collect_base_irf(2,:,i) = transformIRF_used(i)*IRF_extin_bands(i,:,1) ;
    collect_other_irf(i,:,sample) = transformIRF_used(i)*IRF_extin(i,:) ;

end


end

%% plot final figure after loop


figure(6)

subplot(2,1,1)

fill([T_step,fliplr(T_step)],[ collect_base_irf(1,:,1), fliplr(collect_base_irf(2,:,1))], ...
                        grey,'EdgeColor','none', 'FaceAlpha', 0.5), hold on

a = plot(T_step, collect_other_irf(1,:,1),     'k', 'Linewidth', 2);  hold on
b = plot(T_step, collect_other_irf(1,:,2),     'k:', 'Linewidth', 2);  hold on
c = plot(T_step, collect_other_irf(1,:,3),     'k--', 'Linewidth', 2);  hold on
d = plot(T_step, collect_other_irf(1,:,4),     'k-.', 'Linewidth', 2);  hold on
e = plot(T_step, collect_other_irf(1,:,5),     'k', 'Linewidth', 1);  hold on

    plot(T_step,zeros(T_irf,1), 'k')
title('Intervention','FontSize',font_size)
ylabel('USD bn','FontSize',font_size)

subplot(2,1,2)

fill([T_step,fliplr(T_step)],[ collect_base_irf(1,:,2), fliplr(collect_base_irf(2,:,2))], ...
                        grey,'EdgeColor','none', 'FaceAlpha', 0.5), hold on

a = plot(T_step, collect_other_irf(2,:,1),     'k', 'Linewidth', 2);  hold on
b = plot(T_step, collect_other_irf(2,:,2),     'k:', 'Linewidth', 2);  hold on
c = plot(T_step, collect_other_irf(2,:,3),     'k--', 'Linewidth', 2);  hold on
d = plot(T_step, collect_other_irf(2,:,4),     'k-.', 'Linewidth', 2);  hold on
e = plot(T_step, collect_other_irf(2,:,5),     'k', 'Linewidth', 1);  hold on

    plot(T_step,zeros(T_irf,1), 'k')
    
title('USD/JPY','FontSize',font_size)
ylabel('%','FontSize',font_size)

h= legend([a b c d e],'1991-2017','1991-2010','1995-2017','1995-2010','1991-2004') ;
      set(h,'Location','southeast')
    



